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We propose an efficient procedure for determining phase diagrams of systems that are described by 
spin models. It consists of combining cluster algorithms with the method proposed by Sauerwein and 
de Oliveira where the grand canonical potential is obtained directly from the Monte Carlo simulation, 
without the necessity of performing numerical integrations. The cluster algorithm presented in this 
paper eliminates metastability in first order phase transitions allowing us to locate precisely the first- 
order transitions lines. We also produce a different technique for calculating the thermodynamic 
limit of quantities such as the magnetization whose infinite volume limit is not straightforward in first 
order phase transitions. As an application, we study the Andelman model for Langmuir monolayers 
made of chiral molecules that is equivalent to the Blume-Emery-Griffiths spin-1 model. We have 
obtained the phase diagrams in the case where the intermolecular forces favor interactions between 
enantiomers of the same type (homochiral interactions). In particular, we have determined diagrams 
in the surface pressure versus concentration plane which are more relevant from the experimental 
point of view and less usual in numerical studies. 
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I. INTRODUCTION 

The importance of numerical simulations in Physics 
is due to the fact that very few models can be exactly 
solved. In principle one may directly simulate any model 
on a computer. Moreover, the Metropolis [l[ and the 
Glauber [2| algorithms used in Monte Carlo (MC) sim- 
ulations are very general and easy to implement. In 
practice things are not so simple. Near second-order 
phase transitions the configurations generated by these 
algorithms present strong temporal correlations (critical 
slowing down) , which prevent an efficient sampling of the 
configuration space. In addition, hysteresis effects due 
to metastability prevent a precise location of first-order 
transition lines. 

In the last years, several techniques have been pro- 
posed to circumvent these problems such as the reweight- 
ing technique by Berg and Neuhaus @ and the simu- 
lated tempering by Marinari and Parisi A different 
approach is the use of cluster algorithms pioneered by 
Swendsen and Wang [j| and by Wolff @ . Several studies 
have shown the efficiency of cluster algorithms in reduc- 
ing the critical slowing down ;5f. More recently, it has 
been shown that the cluster dynamics may practically 
eliminate metastability in first order phase transitions 
So far this has been achieved only for the Blume-Emery- 
Griffiths (BEG) spin-1 model @, H|, for which a special 
cluster algorithm has been developed [t| [l(| ■ 
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In this paper, we present a simple cluster algorithm 
that eliminates metastability in first-order phase transi- 
tions and combine it to the Sauerwein and de Oliveira 
(SO) method that allows us to obtain the surface pres- 
sure directly from numerical simulations [Tl| . In the orig- 
inal formulation of the SO method, the authors used the 
Metropolis dynamics to generate the system configura- 
tions. However, as explained above, this is not the best 
choice near phase boundaries. We also introduce a simple 
procedure to calculate the order parameters and the con- 
centrations of molecules in the neighborhood of a first- 
order line from numerical simulations. 

As an application of our method, we have determined 
the phase diagrams of the Andelman model for Langmuir 
monolayers made of chiral molecules. More specifically, 
we were interested in surface pressure versus concentra- 
tion phase diagrams, which are more interesting from 
the experimental point of view. The Andelman model 
was first studied by using the mean field approach on a 
bipartite lattice [12j. However, X-ray diffraction exper- 
iments suggest that the condensed phases of Langmuir 
monolayers tend to form triangular structures, and not a 
bipartite lattice as considered by Andelman. In order to 
be more consistent with the physics of Langmuir mono- 
layers, Pelizzola et al [I3| have studied the heterochiral 
case on a two dimensional triangular lattice, using the 
cluster variation method, and have obtained phase dia- 
grams that are qualitatively different from Andelman's. 
We study through MC simulations the remaining ho- 
mochiral case, and we show that, in contrast with the 
heterochiral case, the MC and the mean field methods 
give results that are in good agreement. 

This paper is organized as follows: in section 2 we 
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briefly review the Andelman model, in section 3 we 
present the cluster algorithm and briefly review the 
Sauerwein de Oliveira method, in section 4 we discuss 
the numerical results and conclude in section 5. 



II. CHIRAL LANGMUIR MONOLAYERS AND 
THE BLUME-EMERY-GRIFFITHS MODEL 

Langmuir monolayers are formed by spreading am- 
phiphilic molecules in an air- water interface. Am- 
phiphilic molecules are strongly asymmetric, constituted 
by two parts with opposite features. The first part — 
the head — , is hydrophilic. It is made of polar chemical 
groups and remains on the water. The second part — the 
tail — , is hydrophobic and made of hydrocarbon chains 
which remain in the air. When the tail is strongly hy- 
drophobic, so that the molecules are insoluble in water, 
Langmuir monolayers form a quasi-two dimensional sys- 
tem. This system can be described in terms of the surface 
pressure and temperature and it displays several phases 
with different structural properties (sec, for instance Ref 

A chiral molecule exists in two forms + and — , called 
enantiomers, related by a spatial transformation that in- 
volves a change of parity. An important feature of the 
physics of chiral Langmuir monolayers is the determina- 
tion of the chiral discrimination, which occurs when the 
interaction energy between enantiomers with the same 
chirality is different from the interaction energy of enan- 
tiomers with different chirality. When the intermolecu- 
lar forces favor the attraction between enantiomers of the 
same species, they are denominated homochiral and they 
lead to chiral segregation. On the contrary, if the attrac- 
tion between different enantiomers is favored, they are 
named heterochiral and they lead to a racemic mixture. 

To study the effect of the chirality theoretically, An- 
delman proposed a simple lattice gas model that can be 
described by the Hamiltonian [l3| 



/V -N 



S.J 



(i) 



<i,j> r,s 



where the first sum is over nearest-neighbor pairs, the 
letters i and j denote the sites of a two-dimensional tri- 
angular lattice, the letters r and s denote the enantiomer 
species (r, s = + or — ), e rs are the coupling energies 
(e++ = e__ and e_| = e |_), N r ^ = 0, 1 are the occupa- 
tion numbers at site i, and /i s is the chemical potential 
of the species s. 

This model is equivalent to the Blume-Emery-Griffiths 
(BEG) spin-1 model @, [H, as it can be seen by relating 
the occupation numbers and the spin-1 variables through 
the relations 



N. 
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where a.- L = 0, ±1. Thus, cr, — 1 (—1) represents a + 
(— ) enantiomer and cr,; = a vacancy. In this way, we 



obtain the BEG Hamiltonian 
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The case J > corresponds to the homochiral case (fer- 
romagnetic BEG). When J < 0, we have the heterochi- 
ral one (antiferromagnetic BEG). We will concern our- 
selves with the homochiral case, since it has not been 
studied beyond the mean-field approach. The parame- 
ters J, 4> depend on the interaction energies e++ = e 

and = e |_ between nearest-neighbor enantiomers 

through the formulae 
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The fields H and A are related to the chemical potential 
of the species + and — and they are given by 



H = 



fi + - [i- 



and 



M- 



(6) 
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They are the conjugate parameters of the chiral order 
parameter and the density of enantiomers defined by 



M 



and 



Q 
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(8) 
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where the N± are the total number of enantiomers ± 
and V = L 2 is the number of lattice sites. In particu- 
lar, we are interested in determining the concentration of 
enantiomers + or — (x+ or x_ , respectively) given by 



x± 



(N + ) + (JV_) 
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(10) 



where m = M/V, q — Q/V. 

In this work we shall study first-order transitions be- 
tween concentrated phases where the enantiomers are 
close to each other (phases C+ and C_ rich in enan- 
tiomers of type + and — , respectively) and the so called 
liquid expanded (LE) phases, where there are many va- 
cancies. In the spin-1 language, we shall study transitions 
between ferromagnetic and paramagnetic phases rich in 
zero spins. 
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III. MONTE CARLO METHOD 
A. Cluster algorithm 

In experiments involving chiral Langmuir monolayers 
the non chiral contribution to the interaction energy be- 
tween the enantiomers is usually larger than the chiral 
one. In our simplified model, this corresponds to choos- 
ing the parameter <f> larger than J. In this paper, we will 
consider the ratio <p/J = 3. This choice has also been 
previously made by Andelman [jjj and Pelizzola et al 
[la ]. For <f>/J — 3, Eq. Q can be rewritten, up to a 
constant term, in the following way 



m = -2K £ 5 CTj)CT . + (A - 2Kz) J2<7i-Hj2<Ti 



<»j> 



i=i 



(11) 

where we used the following identity — (er^ix,- + crfaj) — 
2(<x 2 -l)(<xf-l) = -25 ai , aj ,K = pj, A = pA,H = 0H, 
and z is the coordination number. For this Hamiltonian, 
we propose the following cluster algorithm: 

1. Choose randomly a site on the lattice and denote 
CTscod the value of its spin. This is the first spin of 
the cluster (seed). 

2. Choose, with the probability 0.5, one of the two 
other possible spin values that are different from 
<7 scc d- Call this new value <7 new (it will remain fixed 
during the construction of the cluster). For exam- 
ple, if CTgccd is +, cr new can be - or 0. 

3. Activate the links between the seed and its nearest 
neighbors that are equal to cr S eed with probability 
p = 1 — e~ 2A . Each new spin connected to the 
cluster by an activated link is added to the cluster. 
Next, we repeat the activation procedure to all the 
new spins of the cluster. The process stops when 
all nearest neighbors have been tested and no new 
spin is accepted. Now, we attempt to change this 
cluster with spins equal to cr SCG d into a cluster with 

spins (Tnow (see Fig. [TJ for an example of a H > 

transition). 

4. Evaluate the difference SHbuik — T^buik — T^buik, 
where Wbuik is the cluster bulk energy (calculated 
neglecting boundary links) when all spins are equal 
to (Tnow and Hbuik is the cluster bulk energy when 
all spins are equal to cr soc d- If <^buik < 0, we 
change all spins in the cluster to <r ncw with proba- 
bility Pflip(c — > 5") = 1. If <57ibuik > 0, we change 
all spins in the cluster to <7 n cw with probability 
-Pfli P (CT ->■ e) = exp(-/3<57Y b uik)- 

To prove that the algorithm satisfies the detailed bal- 
ance condition we have to consider two types of transi- 
tions: ± <-> =F and ± <-> 0. For the first transition our al- 
gorithm is equivalent to Wolff's @ and for this reason we 
shall concentrate on transitions of the second type. Let 



us consider, to exemplify, the transition 
Eq. (fTI) , we obtain 
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where i ai is the total number of boundary links that 
connect sites with spins a inside the cluster and sites 
with spins 7 outside the cluster. 

The ratio between the transition probability W a —>a 
and the reverse transition probability W&-* a is given by 



W^s _ Wbulk(l -p) £ ++Pfiip(cr -> 5) 



WbuXk(l-pY 00 Pm P (d- -> a) 



(13) 



The bulk term Wbuik is the sum of the probabilities as- 
sociated with all possible ways of activating links to con- 
struct the cluster, the term (1— pY++ is the probability of 
not including in the cluster a nearest neighbor site with 
occupation variable <7 S0C d ■ Analogous comments hold for 
the transition Wa-ta- Clearly, Wbuik = «h>uik because 
for each configuration of activated links in a there is a 
corresponding one in a (see Fig. [1]). Recalling the def- 
inition of Pflip, given in step 4 of the algorithm, we see 
that the ratio of the flipping probabilities is always equal 
to exp(— /3<57ibuik)- Finally, since 1 — p = e~ 2K the right 
hand sides of Eqs. (fl"2"j) and (fi"3"|) are equal and this equal- 
ity implies detailed balance. It is worth mentioning that 
the algorithm proposed here is a particular case of the 
cluster algorithm proposed by Bouabci and Carneiro Q 
and later extended by Rachadi and Benyoussef [l(| for 
other regions of the parameter space. 



B. The Sauerwein and de Oliveira method 

In order to determine the grand canonical potential 
from Monte Carlo simulations, one usually calculates one 
of its derivatives and numerically integrate the results. 
To use this technique one has to know the value of the 
grand canonical potential at a reference point and then 
numerically integrate along a path which connects the 
reference point to the point where one wants to calcu- 
late the grand canonical potential. An alternative is the 
method proposed by Sauerwein and de Oliveira [ll| that 
allows one to directly obtain the grand canonical poten- 
tial from the MC simulation. 

In this method, the largest eigenvalue of the transfer 
matrix is directly evaluated from Monte Carlo simula- 
tions. Since in the thermodynamic limit the grand par- 
tition function is proportional to the largest eigenvalue 
of the transfer matrix, its calculation enables us to de- 
termine all thermodynamic properties, in particular the 
surface pressure that is the negative of the grand canon- 
ical potential. 

In order to explain how to obtain the largest eigen- 
value, let us consider a triangular lattice with V sites 
divided in N successive layers Sk = (&i,k, 02,fe, &L,k) 
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with L spins, V — L x N (All this applies to the triangu- 
lar lattice that we use in our paper.). The Hamiltonian 
may be decomposed in the following way 

N 

H = J2nSk,S k+1 ), (14) 
fc=l 

where due to the periodic boundary conditions Sj\r+i = 
St. The probability P(Si, S 2 , Sm) of a given configu- 
ration of the system is given by 

P{St, S2, — , Sn) — —T(St, S2)T(S2, S 3 )...T(Sn, Si), 

(15) 

where T(Sk, Sk+i) = exp(— /3H(Sk, Sk+i)) is an element 
of the transfer matrix T and 

Z = Tt(T n ), (16) 

is the grand-canonical partition function. By using the 
spectral expansion of the matrix T it is possible to show 
[ill ] that 

<S Su s 2 >=^ <T(St,St)> ■ (17) 

This expression enables us to calculate the largest eigen- 
value Ao of the transfer matrix T in terms of the averages 
< ^Si,S 2 > an d < T(Si,S\) >, where Ss lt s 2 = 1 when 
layers Si and S2 are equal and zero otherwise. We use a 
MC simulation to generate the configurations with which 
we calculate the averages. 

In the specific case of the BEG Hamiltonian in the 
triangular lattice with L x L sites, the transfer matrix T 
of a n-layer is given by 

L 

T(S n ,S n+1 ) = exp{^2[Ka k , n (vk ,n+l "i~ ffe-t-1,71 

fe=l 

+c r fc+i,„+i) + (3^l n (al„ +1 + al +ln + o 2 k+ln+1 ) 
-Aal n +Ha k , n }}. (18) 

The grand canonical potential per site in the lattice gas 
representation (or the free energy in the spin-1 represen- 
tation) is given by 

^ = _JLlnAo = -P, (19) 
where V is the surface pressure. 

IV. NUMERICAL RESULTS 

In this section, we define the following dimensionless 
quantities: 

t = k B T/J, D = A/ J, h = H/J, IT = V/ J, (20) 
where V, the surface pressure, is given by Eq. (|19p . 



As a check on the efficiency of the proposed cluster al- 
gorithm, we show in Fig. [3 the grand canonical potential 
tp versus the chemical potential D for h = and t = 0.8. 
We considered a very low temperature, because in this 
case hysteresis effects are very strong. In Fig. [2] we com- 
pare the performances of the Metropolis and the cluster 
algorithm on a triangular lattice with periodic boundary 
conditions and linear dimension L = 30. To evaluate ip 
and to estimate its statistical error after equilibrating the 
systems we have used 5 x 10 4 Monte Carlo steps divided 
into 1000 independent runs. Note that with the Metropo- 
lis algorithm the system is trapped in metastable states 
and even after millions of MC steps it does not undergo 
a transition to the stable phase. This does not happen 
with the cluster algorithm because the system is able to 
easily pass from one phase to the other. The efficiency 
of the algorithm allows us to determine first-order tran- 
sition lines with high precision and the good quality of 
the data enables us to perform very precise finite size 
analysis. 

In principle it is possible to determine the transition 
point using the free energy. As it can be seen in Fig. [21 
there is a kink in the free energy as a function of D at 
the transition point D* L . One can then perform a finite 
size analysis to obtain D^. However, it is simpler and 
more efficient to analyze the susceptibility whose finite 
size behavior is well known for both first and second-order 
phase transitions. After determining we use the SO 
method at this point to calculate the surface pressure. In 
first-order phase transitions, the surface pressure, which 
is proportional to the negative of the grand canonical 
potential, does not have a finite size behavior as simple 
as the susceptibility. As a function of the system size, 
the surface pressure saturates quickly. Thus, the values 
of the surface pressure that we use in our graphs come 
from the largest lattices that we have simulated. 

The susceptibility is defined as \t = L 2 ((m 2 ) — 
(\m\) 2 )/t, where the magnetization m — J2i a i/V- For 
a fixed system size L, maintaining t and h fixed, and 
increasing D towards the coexistence line, one observes 
a peak in the susceptibility at D* L , as seen in Fig. El 
where the lines were drawn only to guide the eye. In the 
thermodynamic limit, this peak becomes a delta- function 
singularity. According to Refs. [3, [IB], the deviation 
of D* L from its asymptotic value decays as £ -2 , in 
agreement with our results, shown in the inset of Fig. |3]. 
From this law, we have obtained the extrapolated value 
D*^ = 12.0000(1). 

In order to understand this result let us perform an 
exact zero temperature calculation of the transition point 
D* L . At zero temperature the free energy F = U — TS = 
U =< H >, where H is given in Eq. ((3]). The system 
chooses the phase that minimizes the energy U. At T = 
0, for small values of D all spins arc +lif/i>0or— 1 
if h < 0. If D is large enough, the energy is minimized 
when all spins are 0. The transition line is obtained by 
equating the energies in the ferromagnetic (condensed) 
phase with all cr, = + (or — ) and in the paramagnetic 
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(liquid expanded) phase with all cr, = 0. Taking into 
account Eq. ([3]), with 4> = 3J, and the definitions given in 
Eq. ([20| , we calculate the energy U± of the ferromagnetic 
phases and the energy Uq of the paramagnetic phase 



U± = VJ(- 12 =F h + D); U Q = 0. 



(21) 



The equation U± — Uq ^=> D = 12 ± h gives the tran- 
sition lines between the ferromagnetic and paramagnetic 
phases. The equation U+ = U- h = gives the tran- 
sition line between the two ferromagnetic phases (this 
holds for D < 12, for D > 12 the systems is in the 
paramagnetic phase). All these transition lines are rep- 
resented in Fig. [4] that gives the phase diagram of the 
Langmuir monolayer in the plane of the chemical poten- 
tials h x D (Recall that hJ = H = (/i+ — /i-)/2 and 
DJ = A = (— /i+ — /Lt_)/2.). The circles are the results 
of Monte Carlo simulations performed at t — 2.4. The 
error bars are smaller than the circles. It is interesting 
to remark that in the temperature interval relevant for 
Langmuir monolayers the zero temperature calculations 
give practically the same results as the MC simulations 
and the mean-field calculations for the transition lines. In 
Langmuir monolayers language, for h — and low values 
of D (higher chemical potentials), we have the condensed 
phase characterized by a 1 : 1 mixture of the two enan- 
tiomers. For higher values of D a transition from the 
condensed phase, rich in enantiomers ±1, to the phase 
poor in enantiomers, the liquid expanded phase, takes 
place. When the chemical potential of the species are 
different (h ^ 0), we have larger fraction of enantiomers 
+ (— ) whenever h > (h < 0), and in the limit of h >> 
(<< 0) the solution only contains the enantiomer + (— ). 

Another procedure for locating the phase transition 
consists in determining the crossing point of the q versus 
D isotherms for different system sizes. As showed in the 
Ref. [13] , the crossing point is independent of the lattice 
size and properly identifies the transition, as shown in 
Fig. [5l We shall present below an independent derivation 
of this important result based on the work of Borgs and 
Kotecky Q. 

For ft, = all curves of q versus D cross at D* = 
12.0000(1) and q « 2/3 for this value of D. This criterion 
for estimating the value of D for which the phase transi- 
tion takes place agrees very well with the finite size anal- 
ysis of the susceptibility \t that we have discussed above. 
For h ^ 0, two phases coexist at the point D* h which now 
depends on h and all isotherms cross at q » 0.5. We 
remark that if single flip algorithms are used to gener- 
ate the dynamics, one will not be able to determine the 
crossing of the curves due to hysteresis effects. 

More relevant from the point of view of Langmuir 
monolayers, and other physical systems involving mix- 
ture of molecules, are the surface-pressure versus con- 
centration diagrams. But before discussing this phase 
diagram we will describe our procedure to fit the curves 
in Fig. O and to obtain the V — > oo limit of q and m 
that are used to determine the concentration x + (see Eq. 
(JTUJ) ) . Since the simulated system is finite, the calculated 



quantities will be affected by finite size effects. As men- 
tioned previously, in the last years, the finite size theory 
of first order phase transitions has been studied exten- 
sively for quantities, such as the specific heat and the 
susceptibility. There are fewer studies for the dependence 
on the system size of quantities like the magnetization or 
the concentration of molecules [HI, [l!| . In the following, 
we propose a method to determine the concentrations of 
the phases that coexist directly from the numerical sim- 
ulations. The first step consists in noting that q x D (or 
to x D) isotherms can be fitted by the equation 



ce 



-aSD 



(22) 



where a, b, c and d are fitting parameters and 5D = 
D — . We are going to show below that a depends on 
the system size L and the temperature T. An analogous 
expression can be written down for the order parameter 
to. The expression above was inspired by the work of 
Borgs and Kotecky where it is shown that at low 
temperatures the partition function for two coexisting 
phases can be written as 



Z 



-0M0,h)V 



-pf a (P,h)v 



](1 



-L/L \ 



(23) 



where h is the magnetic field (our system also depends on 
the crystal field D), L is a constant of the order of the 
infinite volume correlation length and /$ is a metastable 
free energy for the phase i (i = 1 or 2). 

In our system, for h — three phases coexist at the 
triple point (D^ = 12). Thus, we expect the sum of 
three exponentials instead of two as in Eq. (|23|) . We 
have assumed that all three exponentials have the same 
weight (we shall use our results to check this point). In 
the neighborhood of the triple point 



e -PVfo + e -fiVf+ 



(24) 



where the fa = fi((3, h, D), i = 0, ±, are respectively the 
metastable free energies of the paramagnetic and ferro- 
magnet phases. Away from the coexistence curve, only 
the fi associated with the correct phase remains and be- 
comes the free energy of the system [Z = exp(—f3Vfi) 
)• 

The parameters m and q are given by 



1 dlogZ 
"fJV 3D : 



1 dlogZ 
~pV dh ■ 



(25) 



For h = 0, /+ = /_ = f±. Taking into account this 
fact and Eqs. (|24[) and (f25|) . we can write the parameter 
q as 



(df /dD)e-P v f° + 2{df±/dD)e-P v f± 
e -pvf a + 2e -/3V/± 



(26) 



At the triple point = f±, where /* = fi(0,h 
0, D = D^), i = 0, ± and the exponentials in Eq. 
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which contain the only dependence on the lattice size, 
cancel out and we obtain 



dfo 



dD 



D=D* 



J dD 



D=D* 



(27) 

This is the reason why the q x D curves for different 
lattice sizes cross at the same point. The crossing point 
provides another method do locate the phase boundaries. 

Our calculations are performed at low temperatures. 
In mean-field, the temperatures are measured in units of 
the coordination number z. For the triangular lattice, 
z — 6. Our MC temperature t = 2.4 is equivalent to a 
t = 0.4 mean-field temperature. We can use the exact 
zero temperature energies given in Eq. (|2f [) to estimate 
the derivatives in Eq. (|2"7| . Recalling that fa = Ui/VJ, 
i = 0, ±, at t = 0, we obtain dfo/dD = and df±/dD — 
1. Thus, at the triple point q* w 2/3 which is the result 
that we obtain in Fig. [5] 

An analogous demonstration holds for the case h ^ 0, 
where Z is the sum of two exponentials as in Eq. (|23p . 
The factor 2 in Eq. (|2"9"|) is replaced by 1 and the crossing 
of the q x D curves occurs at the point D « 0.5. 

In the curves plotted in Fig. D varies in the interval 
[11.994, 12.006] which is very narrow. It is possible in 
this case to expand the fa = fi((3,h — 0,D), i = 0,±, 
around the triple point, 



fi = f* + fi'* 5D + 0((SD) 2 ), 



(28) 



where SD = D - D^, f* = f0, h = 0,D = D* oo ) and 
//* = (dfi/dD)\{D = ££,}, for * = 0,±, 



J, zsz, 

/„'* e"' 9l '/o'* 4£ > + 2 /£ e~ /3y/ ± 4£) 

e -/3y/ '* sd + 2e - f3Vf ± SD : 



(29) 



which has the same form as Eq. ([22]) after we divide the 
numerator and the denominator by exp(— /3Vfg* SD). 

In Fig. [5]the symbols stand for the values of q obtained 
from the numerical simulations and the solid lines are fits 
of the points using Eq. (|22p by minimizing the \ 2 merit 
function (20| . In order to perform the fittings we used 
the Levenberg-Marquardt method that is well described 
in Ref . [2(| , where one can also find the subroutines that 
are necessary to implement the method. These subrou- 
tines return the variances of the fitting parameters and 
the quality of the fitting. A few words about the im- 
plementation of the subroutines is in order. Our fitting 
function Eq. (j2"2")l contains exponentials whose arguments 
may become very large. In order to avoid numerical over- 
flow it is convenient to use the asymptotic values of q 
when \SD\ becomes too large. Define, for example, q = b 
for 5D > 30 and q = c/d for SD < —30. Of course, 
the number 30 is rather arbitrary. Non-linear fittings de- 
pend on a good initial guess of the fitting parameters. 
One may proceed as follows. Note that b — q(D — > oo) 
and c/d — q(D — > — oo). Since in the simulations the 
D interval is finite, instead of taking the \D\ — > oo limit 
we use the values of q in our data set associated with 



the largest and the smallest values of D. Call them q + 
and q-, respectively and put b w q + , c/d g_. Next 
define q* = q(D = D*^) and q x = q(D = D\), where 
D\ < D^, is chosen in the region where the graph qx D 
has already started to curve down. It is simple to solve 
the fitting parameters in terms of these quantities. 



1 



D* OQ -D l 



log 

q- (q 



(<?i - q+){q- - q*) 



(q- - qi)(q* - q+) 



q- - q* 



(30) 



With this choice for the initial parameters, the conver- 
gence of the fitting routine is very fast and the quality 
of the fitting is very good (the factor Q that measures 
the goodness-of-fit [20[ is close to 1). In Tables HI and ITT1 
the errors of the parameters are the square roots of the 
variances (standard deviations) that are returned by the 
fitting routines. 

The fitting parameters for the curves in Fig. [5] are 
given in Table Q] Now we can check the equal weight 
hypothesis for the exponentials. For h = the two con- 
densed phases C± have the same free energy and two 
of the three exponentials are identical, as we discussed 
above. The q x D curves in Fig. [5] are in the vicinity 
of the triple point, so we expect that d w 2. This is the 
result that we obtain (see Tables HI and ITT1 for the magne- 
tization to). 

For h ^ 0, there is the coexistence of two phases (LE 
and C+ or C_). Near the transition we have the sum of 
two exponentials with the same weight, as in Eq. (|23[) . 
We have checked that d 1 near the transition line, as 
it was expected. 

Comparing Eqs. (12"2")) and (|2l)|) we note that the param- 
eter a that appears in the exponent should be propor- 
tional to the system volume. A log(et) x log(L) plot gives 
the straight line log(a) = A+B log(L) with A = -0.87(2) 
and B = 1.990(6) for table [fl and A = -0.88(2) and 
B = 1.993(5) for table ITT1 Thus, as expected, the con- 
stant a scales with the volume of the system. 

The values of q for the condensed and liquid expanded 
phases are calculated by taking the L — > oo limit in Eq. 
(|2"2")l . The condensed phases occur in the region D — 
D^ < and for this reason q — > c/d as L — > oo. The 
liquid expanded phase occurs in the region D — D*^ > 
and q — > b as L — > oo. The curves for to x D are very 
similar to the curves q x D in Fig. [5] and can be also 
be fitted by an expression analogous to Eq. [22] Having 
calculated q and to, we obtain x + through the expression 
x + = (1 + \m\/q)/2. The use of \m\ instead of m is due 
to technical reasons (see sections 2.3.3 and 2.3.4 in Ref. 
[2l|). As a consequence, the magnetization is small but 
not zero when h = 0. This introduces a small distortion 
in the diagram of Fig. [S]ncar x + = 0.5, but symmetry 
arguments guarantee that m = when h = and the 
coexistence curve passes through the point with x+ = 0.5 
(filled circle in Fig. [S]). 

The h > half side of the diagram h x D is mapped 
onto the right hand side of the surface-pressure versus 
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concentration diagram (x + > 0.5) whereas h < corre- 
sponds to the x + < 0.5 concentration range. As men- 
tioned above, h = implies that the fraction of enan- 
tiomcrs + and — are equal and in the coexistence of the 
three phases, one has x+ = 0.5. From the point of view 
of homochiral Langmuir monolayers, the chiral segrega- 
tion takes place, in contrast to the heterochiral case, in 
which one has a racemic mixture. The surface pressure of 
a 1:1 mixture of enantiomers is higher than the pressure 
for pure enantiomers. This feature is verified in experi- 
ments in which the chiral segregation occurs [Til ] . Unfor- 
tunately, to date few experiments have been performed 
covering the whole range of concentrations, usually they 
are restricted to the 1:1 mixture and the pure cases. For 
comparison, we have also plotted in Fig. [6] the results 
obtained from the mean field technique. 

In contrast to the heterochiral case, for which the mean 
field results disagree with those obtained from the clus- 
ter variational method, in the homochiral case the accor- 
dance between mean field and the numerical simulations 
is very good. 

V. CONCLUSIONS 

In this paper, we present an efficient way for deter- 
mining phase diagrams from numerical simulations. To 
illustrate it, we have considered a simple model that de- 
scribes the behavior of homochiral Langmuir monolayers, 
which is equivalent to the BEG model. It is worth men- 
tioning that although we have interpreted the phase di- 
agrams obtained here in terms of Langmuir monolayers, 
similar phase diagrams are obtained when one uses the 
BEG model to describe a mixture of two distinct species 



with vacancies. The use of a cluster algorithm that elimi- 
nates metastability in first order phase transitions allows 
us to precisely locate the first-order transitions lines. To 
determine the surface pressure we used the method pro- 
posed by Sauerwein and de Oliveira in which the surface 
pressure is determined directly from the numerical sim- 
ulations without the necessity of performing numerical 
integrations. The fitting procedure proposed in this pa- 
per to determine the concentrations, based on the work 
of Borgs and Kotecky [3] , is easy to implement and uses 
all information contained in the order parameter curve. 
It seems to improve on the usual finite size analysis for 
the magnetization near first-order transition lines in that 
it does not present "overshooting" effects [HI, HH and 
both m and q present a monotonic behavior as a func- 
tion of L, but this point has to be further investigated by 
increasing the statistics. The elimination of metastabil- 
ity also enables us to use the crossing of the curves qx D 
(or m x D) for different lattice sizes as a criterium for 
locating the phase boundaries. This usually cannot be 
done due to hysteresis effects. Finally, we remark that 
our approach is general and it can be used for any spin 
model. In systems for which a cluster algorithm is not 
available, we can use other techniques, such as the multi- 
canonical approach [l| or the simulated tempering [i[ to 
generate the dynamics. 
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FIG. 1: Example of a possible cluster transition. The heavy 
lines are the activated links. 
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FIG. 2: Grand canonical potential tjj/J versus D across the 
first-order line for L — 30, h = and t — 0.8. The first 
graph refers to the Metropolis algorithm and the second to 
the cluster algorithm. The symbol o ( x) indicates increasing 
(decreasing) D. The symbols are larger than the error bars. 
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FIG. 3: Susceptibility \t versus D for several values of system 
size L, h — and t — 2.4. In the inset, we plotted the value of 
D for which the susceptibility is maximum (D* L ) versus L~ 2 . 



L 


a 


b 


c 


d 


24 


233.2(9) 


0.0173(7) 


1.983(8) 


2.000(8) 


30 


363(1) 


0.0166(4) 


1.986(2) 


2.004(9) 


36 


521(2) 


0.0161(4) 


1.99(1) 


2.01(1) 


42 


708(3) 


0.0160(3) 


1.99(1) 


2.00(1) 


48 


925(4) 


0.0160(3) 


1.99(1) 


2.00(1) 


54 


1168(6) 


0.0160(3) 


1.98(1) 


2.00(1) 


60 


1.45(1) xlO 3 


0.0165(7) 


1.98(1) 


2.00(2) 



TABLE I: Values of the fitting parameters obtained from the 
q x D curves in Fig. [4] (i = 2.4 and h — 0). The numbers 
between brackets are the uncertainties in the last digits. 
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FIG. 4: Phase diagram in the space of the chemical potentials 
h versus D. The symbols C+ and C_ denote the condensed 
phases rich in enantiomers + and — , respectively, and LE 
is the liquid expanded phase. The solid line is the t = 
calculation which practically coincides with the mean field 
result. The circles are from MC simulations. 




FIG. 5: Order parameter q versus D for h = 0, t = 2.4 and 
several system sizes L. In the inset, a collapse of all curves 
by plotting q versus z = (D — D^o) * L 2 . 




FIG. 6: Surface pressure II versus concentration x+ phase 
diagram obtained by numerical simulations (circles) and mean 
field (lines). 



L 


a 


b 


c 


d 


24 


233.4(9) 


0.0047(6) 


1.943(8) 


1.998(8) 


30 


364(1) 


0.0036(2) 


1.943(9) 


1.998(9) 


36 


523(2) 


0.0029(1) 


1.94(1) 


2.00(1) 


42 


711(3) 


0.0025(1) 


1.94(1) 


2.00(1) 


48 


930(4) 


0.0022(1) 


1.94(1) 


1.99(1) 


54 


1175(6) 


0.00200(9) 


1.93(1) 


1.99(1) 


60 


1.45(1) xlO 3 


0.0018(3) 


1.93(2) 


1.99(2) 



TABLE II: Values of the fitting parameters obtained from the 
m x D curves for t = 2.4 and h — 0. 



